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ABSTRACT 

I propose a modification of the spherical infall model for the evolution of density fluc- 
tuations with initially Gaussian probability distribution and scale-free power spectra 
in Einstein-dc Sitter universe as developed by Hoffman & Shaham. I introduce a gen- 
eralized form of the initial density distribution around an overdense region and cut 
it off at half the inter-peak separation accounting in this way for the presence of the 
neighbouring fluctuations. Contrary to the original predictions of Hoffman & Shaham 
the resulting density profiles within virial radii no longer have power-law shape but 
their steepness increases with distance. The profiles of halos of galactic mass are well 
fitted by the universal profile formula of changing slope obtained as a result of A^- 
body simulations by Navarro, Frenk & White. The trend of steeper profiles for smaller 
masses and higher spectral indices is also reproduced. The agreement between the 
model and simulations is better for smaller masses and lower spectral indices which 
suggests that galaxies form mainly by accretion while formation of clusters involves 
merging. 

Key words: methods: analytical - cosmology: theory - galaxies: clustering - galax- 
ies: formation ~ large-scale structure of Universe 



1 INTRODUCTION 



I It is generally believed that the structure in the universe 
■ emerged hierarchically from initially small density fluctua- 
$_( ' tions. Small fluctuations, which at present remain such only 
at very large smoothing scales, can be successfully described 
by the linear theory. On smaller scales however, nonlin- 
ear effects come into play and linear approximation is no 
longer valid. In the weakly nonlinear regime perturbative 
approach proved extremely successful in predicting the sta- 
tistical properties of density fields at scales exceeding 10 
Mpc. In order to predict properties of single objects we must 
however trace the evolution of density up to strongly non- 
linear regime. The price to pay is high: we have to restrict 
the analysis to one object neglecting its interactions with 
the neighbourhood. 

The simplest deterministic model of strongly nonlin- 
ear evolution proposed by Gunn & Gott (1972), called the 
spherical model, described the behaviour of a uniformly 
overdense region in the otherwise unperturbed, expanding 
Universe. It was extended by Gott (1975) and Gunn (1977) 
to apply to the evolution of matter around an already col- 
lapsed density perturbation superposed on a homogeneous 
background. The main prediction of the model (called the 
spherical accretion or the secondary infall model, hereafter 



SIM) was that the matter collapsing around the perturba- 
tion should form a halo with r"^''* density profile. 

It is much more realistic to assume that the progeni- 
tors of structure were not the collapsed perturbations but 
the local maxima (rare events) in the density field which 
had initially Gaussian probability distribution. This was the 
approach of Hoffman & Shaham (1985, hereafter HS) who 
applied SIM to the hierarchical clustering scenario. They 
assumed that the density peak dominates to some extent 
the surroundings causing the collapse of the material that is 
gravitationally bound to it. The initial density profile around 
the peak was approximated by the mean density predicted 
by the initial probability distribution which was assumed to 
be Gaussian. HS considered scale-free initial power spectra 
of fiuctuations in difi'erent cosmologies and found that the 
final profiles of halos steepen for higher spectral indices and 
lower density parameter fio- 

This kind of dependence on cosmological parameters 
was confirmed via the results of completely difi'erent ap- 
proach to studies of structure formation: the A^-body sim- 
ulations. Quinn, Salmon & Zurek (1986), Efstathiou et al. 
(1988) and more recently Crone, Evrard & Richstone (1994) 
among others confirmed the analytical predictions by fitting 
power laws to the density profiles of dark matter halos re- 
sulting from their simulations. 

Recently Navarro, Frenk & White (1997, hereafter 
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NFW) performed a series of high-resolution A^'-body sim- 
ulations for power-law initial spectra and found that the 
density profiles of dark halos in large range of masses can 
be fitted with a simple formula with only one fitting param- 
eter. The density profile was observed to steepen from 
near the centre of the halo to at large distances. This 
confirmed earlier results of these authors obtained for CDM 
cosmologies. Similar shapes of profiles were also observed by 
Cole & Lacey (1996) and Tormen, Bouchet & White (1997) 
in their A^-body simulations. 

Although the overall trend of steeper profiles for higher 
spectral indices and lower Qq predicted by SIM was con- 
firmed, NFW claim that the changing slope of the profile is 
"at odds with the analytic predictions" . In fact this is only 
true in the case of fio = 1 for which SIM indeed predicts a 
power-law profile. However, for large (i.e. close to 1) values 
of fio SIM in the form proposed by HS is least reliable. This 
is mainly due to the fact that SIM describes the evolution 
of a single overdense region while in reality the halos form- 
ing at various locations compete for mass. As a result, the 
mass accreted by an overdense region is not the whole mass 
which is gravitationally bound to it (in the case of flo — 1 
this mass is infinite). Another limitation of SIM comes from 
the statistical approach applied in determining the initial 
conditions: the initial density distribution has a finite co- 
herence scale which also may influence the amount of mass 
accreted onto the peak. 

The paper is organized as follows: after a short pre- 
sentation of SIM in Section 2, in Section 3 and 4 I outline 
proposed modifications to the model and in Section 5 com- 
pare obtained halo density profiles to the results of A*'-body 
simulations. The discussion follows in Section 6. 



2 THE STANDARD SPHERICAL INFALL 
MODEL FOR n = l 

In the following I present the main assumptions and results 
of SIM as applied to the density contrast field with initially 
Gaussian probability distribution. This version of the model 
was first presented for the case of linear fields by HS and 
slightly modified by Lokas (1998) in attempt to take into 
account the weakly nonlinear corrections. In the following 
only the linear approximation will be used to determine ini- 
tial conditions. The remaining assumptions and conventions 
will follow those of Lokas (1998) except for a change in nota- 
tion introduced in order to conform to the widely accepted 
notation for the parameters of the universal profile. 

I will consider the Einstein-de Sitter cosmological model 
with n = 1 and zero cosmological constant. The initial prob- 
ability distribution of the density field will be assumed to be 
Gaussian. The initial density power spectrum will be mod- 
eled by the scale-free form 



P(fe) = Cfc" 



(1) 



with indices n — — 1.5,— 1,-0. 5, which are of biggest cos- 
mological interest. The fields will be smoothed with a filter 
of a Gaussian shape 



WG{kR)=e 



(2) 



For such spectra and smoothing the density field can be 
characterized by the correlation coefficient 



n + 3 3 
2' 



X 



(3) 



where ^R{r) is the (auto)correlation function of the two val- 
ues of the density fields (smoothed at comoving scale R) at 
points separated by distance r, x is the distance in units of 
R: X — r/R and a is the linear rms fiuctuation at scale R 
which in this case is given by 



" (2^)2fln + 3 

where the time dependence is 

m 



(4) 



(5) 



1 + z 

Let us assume that at r = we detect an overdense 
region of density equal to a standard deviations: S{r = 0) = 
aa (for high enough a this is approximately equivalent to 
assuming there is a peak at r = 0). Two-point probability 
distribution function then predicts that at distance r (or x if 
we measure the distance in units of R) the expected density 
contrast is 



{S{x)) — g{x)aa 



(6) 



with q{x) given by equation (|^. 

The dynamical evolution of matter at the distance Xi 
from the peak is determined by the mean cumulative density 
perturbation within Xi which is given by 



Ai(xi) 



S{x)x'^dx. 



(7) 



If we assume that S{x) — {S{x)) the expected cumulative 
density can be approximated by 



(Ai(^i)) 



aah{n)x^ 



■(n+3) 



for <C 1 
for Xi ^ 1. 



(8) 



The values of the numerical factor h{n) can be found in 
Lokas (1998). A similar approximation for large Xi was used 
by HS as the initial condition for their SIM. The approx- 
imation at large Xi in the form given above is accurate to 
within 10% for xi > 5. 

The cumulative density contrast Ai(a;i) describes the 
initial density distribution around the peak. Assuming that 
the cumulative density can be approximated by the mean, 
Ai(a;i) = {Ai(a;i)), SIM can be used to predict the final 
density profile which will form after shells which are bound 
to the peak collapse onto it. The distance to the farthest 
shell which is bound to collapse is given by the condition of 
the vanishing energy 



Ai(xi,o) 



: fir 



1 = 



(9) 



where fli is the density parameter at some initial epoch at 
which we specify the initial density distribution. In the case 
of fi = 1 the scale Xifi of the gravitational infiuence of any 
overdense region is infinite, the region should therefore col- 
lect mass from the whole universe. 

According to SIM after virialization the shell Xi ends up 
at X = fxm where Xm is the maximum radius given by 



Xni — Xi- 



Ai + 1 
Ai 



(10) 



and / is the collapse factor, which in the case of scale free 
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initial conditions for large xi, equation (|8|), can be approx- 
imated by a constant of the order of 1/2 (see Zaroubi & 
Hoffman 1993). The final density profile then follows from 
the conservation of mass 

p{x)x^dx — pi{x{)xlAxi. (11) 

If we approximate the initial density of the shell of radius Xi 

Pi(a;i) = Pb,i[l + 5i(a;i)] (12) 

by the background (critical) density pb,i = Pcrit.i = (1 + 
2:i)^Pcrit,o and expand the right hand side of equation ( [lo[ ) 
in Ai keeping only the linear term, we will end up with the 
power law density profile found by HS 

3/{n+4) 

X~ 



pixl ^ (l + z,f 
n + 4 



Pcrit.O 



aoh 

T 



-3(n + 3)/(n+4) ^^^^ 



The density profiles will be expressed here as above in units 
of the present critical density to ensure direct comparison 
with the results of A'^-body simulations. 

It is interesting to note that what is usually quoted as 
the prediction of SIM for hierarchical clustering scenarios is 
the special case ( |l^ which is valid only for f2 = 1. In the 
case of open universe the slope of the halo steepens gently 
with the distance from the centre starting with the slope 
(^^ near the centre and approaching x~* for shells close to 
the last bound shell. The case of f2 < 1 will be treated in 
the follow-up paper. 



3 HOW TO IMPROVE THE SPHERICAL 
INFALL MODEL? 

The main and most questionable assumptions underlying 
SIM as formulated by HS are of course the spherical sym- 
metry of the problem and the absence of peculiar velocities. 
I will keep the assumptions here and show that even with 
these simplifications the agreement between the model and 
the results of A^'-body simulations can be improved. 

First I propose to relinquish the large Xi approximation 
for the (Ai(a;i)), equation (^) applied by HS in their cal- 
culations. The general formula for the expected cumulative 
density contrast can be found using equations (^), (^) and 

(i) 



6aa 



iFi 



{n + l)xf 

'n + 1 1 
2~' 2 



1 ^ 

'2' 



x^ 



iFi 



Xj 

'T 



(14) 



The main reason for this generalization is that espe- 
cially with the conditions given below the main contribution 
to the final density profile of the halo comes from xi of the or- 
der of few. Although HS used the large Xi approximation as 
the one applicable at large distances from the peak, it is pos- 
sible to interpret large Xi also as small smoothing scales. The 
smoothing procedure of the initial density field that is per- 
formed here in order to determine the location of the over- 
dense regions has no equivalent in real structure formation 
or A'^-body simulations where the only artificial scales are 
those of softening length and the size of the simulation box. 
One may argue that the large distance (or small smoothing 
scale) approximation to {Ai(a;i)) is therefore more realistic. 



Table 1. The values of the scales of influence cqh a-ud x; pp/2 
for peaks of height a = 3 for different spectral indices n. 



Xi.cOH 



p/2 



-1.5 14.3 

-1.0 9.02 

-0.5 6.66 

0.0 5.35 



7.36 
6.45 
5.81 
5.32 




a 



Figure 1. The scale of coherence xi cofl (dashed lines) and half 
the distance between peaks x; pp/2 (solid lines) as functions of the 
height of the peak a for different spectral indices. The four curves 
in each set correspond from top to bottom to n = —1.5, —1, —0.5 
and n = 0. 



However, taking into account the intrinsic role of the initial 
smoothing in determination of the initial conditions for col- 
lapse, the masses of halos etc. it is difficult to accept such 
argument . 

Another limitation of the model presented in Section 2 
comes from the statistical approach to determining the ini- 
tial conditions for collapse. The density profile is evaluated 
with the assumption Ai(xi) — (Ai(a;i)), while this approx- 
imation can only be considered valid for scales up to the 
scale of coherence a;i,coH defined by 



(Ai) = Ea. 



(15) 



The calculation of the dispersion Ea in the case of power- 
law spectra and Gaussian smoothing is presented in Ap- 
pendix A. Equation ( p^ can be then solved numerically 
for a;i,coH. The results for different spectral indices n and 
heights of the peak a are shown in Figure 0. Table ^ lists 
the values of a;i,coH for the peak of height a = 3 that will 
be considered in the calculation of the profiles. 

Similar calculations, but without providing the analyti- 
cal expressions for (Ai) and Ea, were performed by HS and 
also by Peebles (1984) and Ryden (1988) who considered 
peaks instead of overdense regions. HS and Ryden (1988) 
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then consider the fraction of mass subscribed to peaks higher 
than aa in the case of different power spectra 



F(> a) = 



47r 



[n,coH(i^)]^iVpk(i')di^ 



(16) 



where Np]^{i') is the differential number density of peaks 
given by equations (4.3)-(4.10b) of Bardeen et al. (1986, 
hereafter BBKS). Since ri,coH(i^) grows faster with v for 
lower spectral indices, for example F{> 3) is below unity 
for higher spectral indices to a few for lower spectral in- 
dices. This leads HS to assume that the most probable 
peaks around which structure will form are those for which 
F{> a) = 1. However, this way they end up with a rather 
surprising value of a > 4 for the height of the most probable 
peak for n = —2. From the statistics of peaks by BBKS we 
know that Np]^{i') is a strongly decreasing function for large 
u and peaks with a > 4 are very rare. 

I propose here to consider the distance between peaks 
as an additional constraint on the cumulative density distri- 
bution around a peak. It is reasonable to assume that only 
the peaks higher than a should be considered as "dangerous" 
for the structure forming around our chosen peak. Since the 
number density of peaks higher than a for sufficiently high 
values of a (a = 2, 3, 4) is a strongly decreasing function of 
a, we may assume that the peaks nearest to our peak are of 
height close to a. The scales of influence of our peak and the 
neighbouring one will be similar and we may approximate 
this scale as half the distance between peaks, a;i,pp/2. The 
scale 2^i,pp is determined by the number density of peaks 
higher than a: 



Xi.pp — 

where 



[npk(> a) 



-1/3 



R 



n 



pk 



(>a) 



Npi,{u)du. 



(17) 



(18) 



The numerical results for 2;i,pp/2 for different spectral 
indices are shown in Figure 0. The results for a = 3 are 
also listed in Table |l| In the following I will simulate the 
constraints on the cumulative density ( p^ by cutting it off 
at the smaller of the two scales 2;i,pp/2 and a;i,coH. I have 
chosen to consider here peaks of given height a = 3 which are 
high enough to collapse by themselves and frequent enough 
to produce large number of objects. We see that the scale 
x\,pp/2 for a = 3 is always smaller than the corresponding 
a^i,coH except for the highest spectral index considered, n — 
0, where the two scales are almost equal. This motivates the 
introduction of the cut-off at a:i,pp/2 and not Xi,coH. 

The remaining issue is the shape of the cut-off function. 
I will adopt a sharp cut-off such that the corrected cumula- 
tive density will be 

(Ai(^i)) 



Ai,cut(a;i) = 



(19) 



1 -\- e(a:i-a:i,pp/2)/™ 

with {Ai(a;i)) given by equation (p^. Unfortunately, in con- 
trast to the problem of the scale of cut-off, we cannot pre- 
dict the width of the filter, w from the model. It will be 
treated as an unknown parameter which has to be specified 
by matching the results from SIM and A^'-body simulations 
(see Section 5). 

Summarizing, the final density profile of the virialized 
halo is given by 



p{x) 



= {l + aag)(l + zi 



Pcrit.O 

with 

_ xi/[Ai,cut(a;i) + l] 



3 rxiy dxj 



(20) 



(21) 



Ai,cut(a;i) 

Formula ( po| ) gives a complicated but analytical expression 
for the density profile as a function of the initial radius Xi 
which is related to the final radius x by equation (|2l|). Since 
the initial density distribution in its generalized form ( p^ 
is not scale free the collapse factor / in equation ( pl| ) can 
no longer be assumed to be the same constant for all shells 
(for the results with / = 1/2 see Lokas 1999). Detailed cal- 
culation of this factor is presented in the following section. 



4 THE COLLAPSE FACTOR 

The purpose of this section if to find a correction to the fidu- 
cial density profile defined as the density distribution with 
all shells stopping at their maximal radii. This corresponds 
to putting / = 1 in equations (po|)-(|2]]). In reality after 
reaching the maximum radius a given shell will start oscil- 
lating and it will contribute to the actual mass contained 
inside inner shells which are expected to contract adiabati- 
cally (see the discussion in Zaroubi & Hoffman 1993) by a 
factor 

M{xi) 



f{xO = 



(22) 



M{xi) -|-madd(a;i) 
where M{xi) = M{x) is the mass contained inside the shell 
of radius Xi or x and madd(a^i) is the mass contributed by 
the outer shells. 

If the shape of the fiducial profile is a power-law given 
by ( |l3| ) the collapse factor is easily determined to be 

1 



/ = 



(23) 



1 (3 - 7) «2-7P(l/u)dM 

where 7 = 3(n -I- 3)/(n -I- 4) , U is the size of the system in 
units of the radius of the considered shell, P is the fraction of 
time a particle belonging to shell of radius x' spends within 
radius x 

dx 



P{x,x') 



(24) 



In the case of scale-free profile P can be expressed as a 
function of a single variable it — x/x' . The normalization of 
P is provided by the obvious condition P{x,x' = a;) = 1. 
The velocity v is calculated by integrating the equation of 
motion of the sheU d^r/dt^ = -GM/r^. 

In the scale-free case it is possible to obtain analytical 
solutions for P 



P{u) = 



7r[7/(4~27)]M 
2V^r[l/(2-7)] 



(25) 



2F1 



,2-7 



2' 2 



for 7 < 2 and 

2r[l/(7-2)]u^/2 



P{u) 



V¥7r[7/(27 - 4)] 
'1 7 , 



(26) 



2F1 



2' 27 - 



27- 
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for 7 > 2. These results provide generalizations of the lim- 
iting cases considered by Zaroubi & Hoffman (1993). The 
collapse factor can then be evaluated numerically for any 7 
and U . 

When the fiducial density profile is scale-dependent as 
in the case of interest here the calculation of P and madd 
is more complicated and has to be done numerically. The 
velocity in equation (p3) is 



V oc [E{x') - E{x)] 



1/2 



where 



£(,)0C / ^dX: 



XiAf da; 
H- Ai da^ 



da;i 



(27) 



(28) 



(29) 



where we used the expression for mass within x\ 
M{x) = M{xi) = ^p,,,,^o{xiRf[l + Ai(xi)] 

and equation (^l|) with f = 1. The cumulative density distri- 
bution Ai is taken in generalized form with cut-off, equation 

Due to the complicated form of the integrand in the 
expression (^sj) for E the integration cannot be performed 
analytically, but the integrand can be exactly approximated 
by a polynomial in Xi making the integration possible. Af- 
ter changing variables from x to xi in equation ( p^ ) the 
value of P can be calculated numerically for any {xi,x[) 
pair. The added mass in equation ( p2| ) can then be obtained 
by summing up the input from j intervals between a given 
shell and the shell that is presently turning around Xi^ta (the 
turn-around radius is obtained in a similar way as the virial 
radius Xi^^, from tu = tcoii/2, see the next Section) 



47rpcrit,()-R^(l + Zi)^ 



E 



dyi. 



(30) 



The result for / is sensitive to the number of steps taken 
in the integration but usually j = 10 is enough to obtain / 
with two digit accuracy for Xi > 1. 

After calculating the factor / for a number of Xi values 
we can approximate the function fjxj) by a polynomial in 
Xi and introduce it into equation (|2^). From equation ( |20[ ) 
we will then obtain the final density profile of a halo. 



5 COMPARISON WITH THE UNIVERSAL 
PROFILE 

NFW performed a series of A^-body simulations of structure 
formation in fiat and open universe for different scale-free 
power spectra of the form (|l|). They concluded that the den- 
sity profiles of dark matter halos are well fitted in all cases 
by a formula 



(31) 



pjx) _ (5char 

Pcrit.o (r/rs)(H- r/rs)2 

with a single fitting parameter 5char, the characteristic den- 
sity. The so-called scale radius is defined by 

rs = - (32) 



where Vv is the virial radius i.e. the distance from the centre 
of the halo within which the mean density is v times the crit- 
ical density, c in equation (^^ is the so-called concentration 
parameter, which is related to the characteristic density by 



^char — 



3[ln(l c) - c/(l + c)] ■ 



(33) 



NFW assume v = 200 for all considered cosmological 
models, which is not strictly correct. For = 1 the exact 
prediction of the spherical model for the ratio of the density 
to the critical density for objects collapsing at the present 
epoch is « = IStt'^ « 178 (see e.g. Lacey & Cole 1993) so 
the natural choice would seem to adopt the value v = 200 in 
order to keep the assumptions as close as possible to those 
of NFW. Ifowever, this value is inconsistent with the more 
advanced SIM considered here: the value of v can be calcu- 
lated and it proves to be much lower than 200: it is of the 
order of 30 for the small masses and reaches 200 only for 
the largest halos. Fitting the NFW formula with v = 200 
to the results of SIM produces a significant offset between 
the fit and the SIM results, one therefore has to adopt the 
true value of v. Then the problem of comparison with the 
simulations arises: what is the relation between masses of 
halos determined from simulations with different v7 Fortu- 
nately, at the virial radii the density profiles of halos in the 
simulations have slope —3 so the mass grows very slowly 
(logarythmically) with r and the masses determined with 
V = 200 and v — 30 should be almost the same. 

Since the measurements of the halo density profiles in 
the Ai'-body simulations of NFW are all done at the final 
epoch which is identified with the present, in the following 
I will assume that for all the halos the collapse time is the 
present epoch. According to the spherical model the collapse 
time of the shell Xi in the = 1 universe is (Padmanabhan 
1993) 



trnU — 



7r(l + Ai) 



iyo(l + 2i)3/2Af 



3/2 ■ 



(34) 



The present age of the universe is tu = 2/{3Ho) and solv- 
ing equation fcoii = tu for Ai gives us the values of the 
cumulative density contrast as a function of Zi, which will 
be denoted by Ai,„. This density contrast corresponds to 
the presently collapsing initial shell a;i,„ which can then be 
found by numerically solving equation 



Ai,cut(a;i) = Ai,!,, 



(35) 



and the corresponding final virial radius x^ is obtained from 
equation (^l|). In order to solve this equation we have to 
specify the initial conditions: the height of the peak a and 
the rms fiuctuation of the density field a. I have already 
chosen a = 3 and will assume that for all halos the starting 
point is such that aa = 0.1, small enough for the linear 
theory to be valid. Final results are not very sensitive to 
this particular choice. 

We also need to make connection to the present magni- 
tude of fiuctuations so I will adopt the standard normaliza- 
tion of the density field such that ag — 1 (rms fluctuation in 
spheres of radius 8h~^ Mpc is one). This normalization can 
be also expressed in terms of the present nonlinear mass M* 
used by NFW. This mass is defined by the condition of the 
rms fiuctuation at this mass scale being equal to the present 



© 0000 RAS, MNRAS 000, 000-000 



6 Ewa L. Lokas 



critical threshold for collapse 5crit- Given the dependence of 
(J on the smoothing scale this yields 

6/(n+3) 



M. = M{R) 



(Jth{R) 



(36) 



where R is the smoothing scale used for normalization, 
subscript TH refers to the top-hat smoothing, equation 



(|A7|), and M(R) = iirR'^pb/S. For O = 1 the threshold is 

Mpc we 



Scrit ~ 1.69 so with our normalization at R = 8h 
have 



5.96 (1.69)-'^''<"+=*> X 10^* h-^Me 



(37) 



Once the normalization is set and the conditions a — 3 
and aa — 0.1 are adopted choosing the initial redshift Zi for 
a given spectral index n gives us the comoving smoothing 
scale R with which the overdense regions are identified. The 
mass of the halo within the virial radius a;„ can then also be 
determined. According to the results of the previous section 
the total mass inside the virial radius is 



M = M(a;i.„) +madd(3;i,„). 



(38) 



One of the main results of NFW was the dependence 
of the shape of the density profiles of halos on their mass: 
their profiles were steeper for lower masses. On the other 
hand, the standard prediction of SIM, equation (p^, gives 
the same profile independently of mass. However, with the 
improvements introduced in Section 3 and 4 it is possible to 
reproduce the dependence of the profiles on mass. 

It has been suggested in the literature (e.g. Katz, Quinn 
& Gelb 1993) that A''-body simulations indicate that if the 
density field is smoothed with one scale R lower peaks 
end up as galaxies and higher ones as clusters. This, how- 
ever, would violate the hierarchical way of structure forma- 
tion since higher peaks collapse earlier. Another argument 
against such assumption comes from the calculations based 
on the improved SIM for n — —1: the reasonable range of 
peak heights a between 2 and 4, which are most likely to pro- 
duce halos, leads for a given smoothing scale to the range of 
masses spanning only one order of magnitude, while in TV- 
body simulations NFW observe halos with masses spanning 
few orders of magnitude. This suggests that the dependence 
on mass should rather be related to the initial smoothing 
scale. 

Calculations show that the mass of a halo increases with 
the smoothing scale up to the scale for which the solution 
of equation (M) yields rather small values of Xi^^ (of the 
order of unity). In these cases only the nearest neighbour- 
hood of the peak determined with the smoothing scale R 
managed to collapse by the present. Such peaks cannot be 
considered as collapsed structures because they accreted up 
till now only a small fraction of the mass that is bound to 
them (i.e. contained inside the cut-off scale). Therefore in 
the following I will deal only with objects that have mass 
increasing with the smoothing scale. For clarity hereafter I 
will consider different smoothing scales for only one height 
of the peak a = 3. 

NFW suggest that the dependence of the shape of the 
halo on its mass is related to the formation time of the halo: 
smaller halos that formed earlier have steeper profiles. Al- 
though I have assumed the collapse time to be the present 
epoch for all halos, this faster collapse of smaller halos can 
nevertheless be observed: smaller masses are obtained from 




Figure 2. The collapse factor / as a function of the initial radius 
of a shell (in units of the virial radius) for spectral index n = — 1. 
The upper and lower curves correspond respectively to the small 
and large mass. 



smaller smoothing scales which correspond to higher initial 
redshifts Zi and from equation (^^ we have that constant 
collapse time leads to lower Ai and therefore higher a;;^^. 
Smaller halos form earlier in the sense that their initial red- 
shifts are higher and by the present epoch their virial radii 
reach the cut-off scale. 

Given the initial conditions as specified above we may 
now proceed to the direct comparison of the profiles ob- 
tained from SIM and the simulations. An important step 
in determination of SIM profiles is the calculation of the 
collapse factor (see Section 4). Figure ^ shows the depen- 
dence of this factor on the initial radius of the shell for two 
different masses in the case of n — —1. The upper curve 
in the figure corresponds to the smallest mass available for 
n = -l (with Zi = 1500), i.e. 0.0017M. = 2.1 x 10"/i"^Mq. 
The lower curve is for the largest mass (with Zi — 43), 
3.3M, = 4.1 X lO"/i"^M0. Note that the initial radius was 
expressed in units of the initial virial radius Xi^v, which is 
much larger for the smaller mass. The squares marked in the 
plot show the points corresponding to O.OlXv. Since the den- 
sity profiles will be considered in the range O.Olx^ < x < x^, 
only the collapse factor values to the right of the marked 
points will be taken into account in the calculations. 

The behaviour of the collapse factor confirms the in- 
tuition one may obtain from the scale-free case, equations 
(p3[)-(p^): the collapse factor is smaller for higher effective 
index of the cumulative density distribution (see the discus- 
sion in Zaroubi & Hoffman 1993). Near the center of the 
halo, where the initial density distribution is flat, the outer 
shells contribute significantly to the inner mass and the con- 
traction is stronger. In more steeper parts of the density dis- 
tribution the outer shells contribute less and contraction is 
weaker. The dependence on mass comes from the fact that 
small mass halos have large virial radii Xi^v reaching steeper 
parts of the distribution and even the cut-off scale, while in 
large mass halos the final density distribution comes from 
shells that initially were quite close to the peak. 
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Figure 3. The density profiles of dark matter halos of mass of the order of 0.02Af» for different spectral indices n in the range O.OlXi, < 
X < Xv The solid lines show the predictions of SIM. The long-dashed ones give the NFW results with their fitted concentrations, while 
the short-dashed curves present the NFW formula (|39|) with concentrations fitted to SIM results. 



Table 2. The parameters characterizing the small mass halos whose density profiles are presented in Figure 



n 






R[h-^ Mpc] 






rv[h~^ Mpc] 


A/[lOi2h-iM0] 


M/M, 


CNFWl 


CNFW2 


-1.5 


350 


0.00805 


0.142 


7.54 


879 


0.356 


1.55 


0.0211 


31.0 


33.5 


-1.0 


600 


0.00469 


0.188 


6.76 


1348 


0.422 


2.58 


0.0208 


59.2 


46.9 


-0.5 


1000 


0.00281 


0.226 


6.21 


2065 


0.466 


3.45 


0.0203 


131 


65.7 


0.0 


1500 


0.00187 


0.268 


5.73 


2861 


0.511 


4.53 


0.0212 


306 


84.4 



Table 3. The parameters characterizing the large mass halos whose density profiles are presented in Figure Ul 



n 




Ai,„ 


R[h-^ Mpc] 






r„[/i-l Mpc] 


M[10^^h-^MQ\ 


M/M» 


CNPWl 


CNFW2 


-1.5 


150 


0.0188 


0.438 


6.50 


301 


0.872 


3.18 


0.432 


15.6 


14.0 


-1.0 


180 


0.0157 


0.625 


5.45 


298 


1.03 


5.55 


0.447 


21.6 


15.1 


-0.5 


200 


0.0141 


0.816 


4.67 


280 


1.14 


7.83 


0.461 


34.9 


16.2 


0.0 


200 


0.0141 


1.02 


4.01 


240 


1.22 


9.83 


0.470 


62.0 


16.8 
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Figure 4. The same as Figure ^ but for masses of the order of 0.5Af«. 



In the comparison between the density profiles pre- 
dicted by SIM to the results of the A'^-body simulations of 
NFW I will use the NFW formula @ in the form 



p{x) vc^ 

Pcrit.o 3[ln(l + c) — c/(l + c)]{cx/xy){l + cxjx^Y 



(39) 



where definition of characteristic density ( |33| ) was used and 
the distances were expressed in units of the smoothing radius 
and denoted by x in order to provide direct correspondence 
with the predictions of SIM. NFW claim that their fitting 
formula is valid in the range 0.01a;„ < a:: < a;„. I therefore 
fit the formula (M) to the points (x, p/pcrit,o) obtained from 
equations (|2o|)-(p1]), spaced uniformly in log {x/xv) in this 
range. For a given spectral index n the single fitting param- 
eter is the concentration c. 

NFW consider density profiles of halos of mass roughly 
in the range 0.01 < M/M, < 10 so the initial redshifts 
will be chosen here so as to obtain similar range of masses. 
Examples of profiles for halos of different mass are shown 



in Figures ^ and ^. Figure ^ shows the density profiles for 
objects of mass of the order of 0.02Mt, which corresponds 
to a big galaxy mass scale, while Figure ^ shows the profiles 
for halos of mass of the order of 0.5M, , which is closer to 
the mass of a galaxy cluster. Such choice of the mass values 
is motivated by the range of masses obtained for different 
spectral indices (see Figure ^ and the following discussion). 
The initial redshifts needed to obtain halos of those masses 
and all the following parameters are given in Tables |2| and 
^ for the small and large masses respectively. 

In agreement with the above discussion of the depen- 
dence of mass on the smoothing scale, smaller masses re- 
quire smaller smoothing scales and the condition aa = 0.1 
translates into their larger initial redshifts Zi. In the case of 
small masses the initial virial radii of the order of 

the cut-off scale Xi^^^/2 which means that the cut-off really 
infiuences the formation of small halos, i.e. in the absence of 
the cut-off their virial radii would be much larger. In the case 
of the large mass halos the situation is different: their virial 
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Figure 5. The dependence of the concentration parameter c on the mass. The dashed Uncs show the results of A^-body simulations 
of NFW, and the solid curves represent the results of SIM. In each set the curves from top to bottom correspond to n = 0,-0.5,-1 
and —1.5 respectively. In the left panel the SIM results were obtained with constant w = 0.6 while in the right panel the values were 
w = 0.14, 0.31, 0.6, 0.95 respectively for n = 0, -0.5, -1 and -1.5. 



radii gnificantly smaller than the cut-ofT scale: this 

scale could not be reached by the present epoch. The proper 
virial radii r„ = x^R/il + Zi) are few times larger for larger 
masses and in both cases agree well with the observational 
estimates of the sizes of the halos of galaxies and clusters. 

The last two columns of Tables ^ and ^ give the concen- 
trations CNFWi and cnfw2 of the NFW formula. The values 
of cnpwi were calculated from the model based on merging 
formalism provided by NFW that describes best the results 
of their A''-body simulations, while cnfw2 are the values of 
concentration obtained by fitting formula (^) to the results 
of SIM. Following the arguments presented above the fits 
were done for the value of v calculated from the SIM. The 
width of the filter w was chosen so that the concentrations 
from NFW and SIM match for the smallest halos in the 
n = — 1 case. 

Figure ^ shows how the concentration c depends on the 
mass of the halo and the spectral index in the larger range of 
masses. The plots cover the largest range of mass available 
from SIM for different spectral indices with the initial con- 
ditions I adopted. The smallest masses were obtained for 
Zi — 1500 and the largest correspond to the point where 
the mass starts to decrease with Zi (the initial virial radii 
are then rather small, significantly smaller than the cut-off 
scale). The values of concentrations from the simulations 
were calculated in the same range of mass as the one ob- 
tained from SIM. 

The solid lines give the concentrations obtained from 
fitting SIM results and the dashed lines are the A'^-body re- 
sults of NFW (one has to remember that those curves are 
fits to points displaying some scatter that were themselves 
obtained by fitting the concentrations to the highly irregu- 
lar density profiles of halos in the simulations). In each set 
the curves from top to bottom correspond to n = 0, —0.5, —1 
and —1.5 respectively. The results for higher spectral indices 



are shown in the smaller range of mass because obtaining ha- 
los of smaller mass for example in the case of n = would 
require initial redshift of Zi > 1500 i.e. initial time earlier 
than the recombination epoch. Similar behaviour - more 
massive and more slowly forming halos for larger spectral 
indices - is observed in the A^'-body simulations of NFW. 
The two vertical dotted lines in the left panel indicate the 
mass scales chosen for detailed comparison shown in Fig- 
ures and Tables 

The left panel of Figure |^ shows the SIM results for the 
constant width of the cut-off filter w = 0.6 which matches 
the concentrations for smallest masses in the case of n = —1. 
There is however no reason why the same width of filter 
should apply for different spectral indices so we can try to 
adopt different widths for different spectral indices. Results 
of matching the SIM and NFW concentrations for the small- 
est masses in all cases are shown in the right panel of Fig- 
ure ^ The adopted widths were w = 0.14, 0.31, 0.6, 0.95 re- 
spectively for n = 0, —0.5, —1 and —1.5. Although we have 
no way of estimating the exact shape of the cut-off function 
from theory, one can expect that the cut-off will be sharper 
for peaks with sharper initial density distribution i.e. for 
higher spectral indices. This guess is in agreement with the 
dependence of the best matching values of w on the spectral 
index. 

The predictive power of the improved SIM suffers also 
from the fact that the discrepancy between the SIM and 
NFW concentrations is much more sensitive to the unknown 
parameter w than to the cut-off scale, which we have esti- 
mated from the peaks formalism. Calculations of the profiles 
with different cut-off scales show that e.g. decreasing the cut- 
off scale gives smaller mass and steeper profile but because 
the mass is strongly affected the corresponding concentra- 
tion from the simulations is also decreased and the factor by 
which both concentrations differ remains roughly the same. 
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On the other hand, decreasing the width of the cut-off filter 
does not significantly affect the mass while the concentration 
of SIM profiles is increased. 

Leaving aside the exact determination of the width of 
the filter, Figure ^proves that SIM predicts the same general 
trend in the dependence of the shape of dark halo profiles 
on their mass: the smaller the halo mass the steeper the 
profile. The dependence of the shape of the halo on the ini- 
tial power spectrum is also reproduced in the wide range of 
masses: concentrations are larger for higher spectral indices 
up to a mass of the order of the nonlinear mass M« where 
SIM approximation seems to break down. We can also con- 
clude from Figure ^ that the agreement between SIM and 
the simulations is better for lower spectral indices. One may 
interpret this result as an indication that for spectra with 
more large scale power and long-distance correlations (low 
n) accretion is a realistic mechanism of growth of fluctua- 
tions. On the other hand, in the case of higher spectral in- 
dices isolated halos tend to appear and even relatively small 
halos form by merging of yet smaller objects. 

Another argument for the reliability of SIM comes 
from the possible (and rather straighforward) application of 
the results presented here to more realistic scale-dependent 
power spectra. Since smaller masses correspond to smaller 
smoothing scales where the effective spectral index of a real- 
istic spectrum is lower, the dependence of the concentration 
on mass according to the SIM should be flatter. This effect 
is also well visible in the simulations of NFW: for the CDM 
spectrum the trend of smaller concentrations with growing 
mass is preserved but the dependence on mass is rather weak 
(especially when we take into account large scatter in the re- 
sults). 



6 DISCUSSION 

The improved SIM presented in this paper provides sim- 
ple understanding of the dependence of the shape of the 
halo on its mass: smaller halos start forming earlier and by 
the present epoch their virial radii reach the cut-off scale 
that accounts for the presence of the neighbouring fluctu- 
ations; more massive halos form later and their virial radii 
are not affected by the cut-off scale, their virialized regions 
contain only the material that initially was quite close to 
the peak identified with the smoothing scale corresponding 
to the mass. 

Although the concentrations of profiles predicted by 
SIM depend on the exact shape of the cut-off function, this 
study leads to rather firm conclusion that the agreement be- 
tween the A'^-body and SIM profiles is the better the smaller 
the mass of the halo and the lower the spectral index of 
the initial power spectrum of fluctuations. If the proflles ob- 
tained from the simulations indeed reflect reality this may 
indicate that galactic halos form mostly by accretion, while 
for cluster size objects merging must be taken into account. 
As suggested by Syer & White (1998) the universal pro- 
file can be formed by sufficiently dense satellites reaching 
the centre of a halo intact and forming a cusp. Since the at- 
tempts to obtain the dependence of the profiles on mass only 
from formalisms describing the merging of halos (Nusser & 
Sheth 1999, Avila-Reese, Firmani & Hernandez 1998) have 
not been fully successful, it seems that the best descrip- 



tion of halo proflles should be provided by a model dealing 
with a combination of accretion and merging (for the dis- 
cussion on the distinction between the two phenomena see 
Salvador- Sole, Solanes & Manrique 1998). This work pro- 
vides the understanding of the origin of universal proflle in 
the case of small masses for which SIM can be considered a 
valid approximation. 

One of the important shortcomings of SIM, that was 
not mentioned here, is the artificial combination of the lin- 
ear and strongly nonlinear regime without taking into ac- 
count the weakly nonlinear phase that may affect the initial 
density distribution before the strongly nonlinear evolution 
takes over. As discussed by Lokas (1998) these corrections 
are of similar importance as the ones introduced by using 
the formalism of BBKS to describe peaks in the density field 
instead of overdense regions. Both effects tend to steepen the 
initial density profiles but are difficult to model analytically 
(for the corrections for peaks see e.g. Ryden 1988). With 
the modifications of SIM introduced here we are interested 
mostly in regions not very distant from the centre of the 
forming halo. In these regions weakly nonlinear corrections 
to the expected overdensity (5) are known only numerically 
and it is difficult to obtain the cumulative density { Ai) . Even 
if we approximate it by some analytical expression we can- 
not proceed because the formula for the final profile is so 
complicated that any perturbative treatment is impossible. 
Qualitatively one may expect that the final proflle will be 
somewhat steeper but since the value of Ai,„ will not be 
changed the solution for the virial radius Xi^v will be lower. 
It follows that also the halo mass will be decreased but it 
will have to be compared to a less massive and therefore 
steeper halo from the simulations, so it is difficult to pre- 
dict whether the agreement between the two results would 
be improved. It should be added that this picture of weakly 
nonlinear corrections does not take into account the paral- 
lel evolution of the rms fiuctuation a itself. As discussed by 
Lokas et al. (1996) its value in the weakly nonlinear regime 
may be close to linear in the case of n = — 1 but may dif- 
fer significantly from the linear prediction for other spectral 
indices. 

No solution to the problem of dark halo formation can- 
not be considered valid without a detailed agreement be- 
tween its predictions and observations. To date several such 
comparisons were performed, in most cases in the form of fit- 
ting the NFW profile to the observed profiles of galaxies and 
clusters that are believed to be dominated by dark matter or 
provide some indication on how light traces mass. Carlberg 
et al. (1997) find that universal profile of NFW provides a 
very good fit to the density profiles of clusters in the CNOC 
survey, while Adami et al. (1998) find that galaxy distribu- 
tion in clusters in the ENACS sample displays a core rather 
than a cusp in the central regions, but the mass distribution 
(Adami, private communication) is somewhat steeper. These 
results suggest that the universal profile indeed provides a 
good description of the density distribution in clusters. In 
the case of galaxies the situation is less satisfying. Kravtsov 
et al. (1998) analysed density proflles of dwarf and low sur- 
face brightness galaxies and found that they are much better 
fltted by a so-called Burkert profile (Burkert 1995) with a 
flatter cusp than in the NFW formula. They also performed 
a series of A'^-body simulations based on a different code than 
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that used by NFW and found similar shapes of galaxies as 
those observed. 



(5) = aag 



(A2) 
(A3) 
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APPENDIX A: THE SCALE OF COHERENCE 

The calculation of the scale of coherence a;i,coH of the cumu- 
lative density contrast Ai can be done in a way analogous to 
the calculation of the scale of coherence Xcoh of the density 
contrast S itself (Lokas 1998). The condition for Xcoh was 



{5) = {{S - (3))^^' 
where 



(Al) 



The results were obtained from the two-point Gaussian 
probability distribution function for fields 5 and 7, from 
which a conditional probability distribution for S followed 
given that at distance r from where 5 is measured there is 
an overdense region where the density is 7 = aa. 

Here we may proceed in a similar fashion and consider 
a two-point probability distribution for variables Ai given 
by equation (Q) and 7. If the condition for 7 is the same we 
have (see also HS, Peebles 1984) 



(Ai) = aEg 
((Ai-(Ai))^) 



E^(l 



(A4) 
(A5) 



The symbols E and Ea refer respectively to the uncon- 
strained and constrained dispersions of the Ai field. The 
new correlation coefficient g' — f'/(aE) is the normalized 
correlation function f'(n) = (Ai(ri)7(0)) . 



The definition of Ai, equation (M), can be rewritten as 
1 



Ai(ri) = 



(27r)» 



5{k)WTu{kri)d^k 



where W^TH(fcri) is the top-hat window function 
3[sin(fcr-i) — fcricos(fcri)] 



W^TH(fcri) 



(fcri)3 



Expression (A6) leads to 



(Ai(ri)) = 



(27r)3cr 



Pfl(fc)W^TH(fcri)d^fc 



(A6) 



(A7) 



(A8) 



where PR{k) is the power spectrum initially smoothed with 
the Gaussian filter 



Pfl(fc) = P{k)e-' ^ . 

The expected cumulative density is therefore 
6Ca 



(Ai) = 
where 



-/i(xi) 



(A9) 



(AlO) 



(sin k ~ k cos k) e 



dk 



i.r-r(^) 



n + 1 3 



(All) 



n + 1 1 



Together with the expression for a, equation (j^, this leads 
to the expression (^) for the expected cumulative density. 
The unconstrained variance of Ai is 



(27r)» 



1 Q/^ 

PR{k)W^n{kn)d\ = 



/2(xi)(A12) 
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where 

hixi) = I fc"-*(smfc- fccosfc)^ e-'^'/'J^d/fc (A13) 

cf + 2 (n-\ \ 2 







i.r'r(^) 



2 V 2 ' 2' 



n - 2 ^ /n-3 1 2\ a^f 1 



n-3 V 2 '2' V 2 n-3 
and the correlation coefficient is given by 

^ -r[(n + 3)/2]^' /2- 

Given the analytical expressions for (Ai), E and (J 
equation ( |l5| ) has to be solved numerically to determine 
a^i = 2;i,coH- The results for different power spectra are 
shown in Figure |^. 
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